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Abstract 

Traveling wavetrains in generalized two-species predator-prey models and two- 
component reaction-diffusion equations are considered. The stability of the fixed points 
of the traveling wave ODEs (in the usual "spatial" variable) is considered. For general 
functional forms of the nonlinear prey birthrate/prey deathrate or reaction terms, a 
Hopf bifurcation is shown to occur at two different critical values of the traveling 
wave speed. The post-bifurcation dynamics is investigated for five different functional 
forms of the nonlinearities. In cases where the bifurcation is supercritical, the post- 
bifurcation behaviour yields stable periodic orbits of the traveling-wave ODEs in the 
spatial variable. These correspond to stable periodic wavetrains of the full PDEs. 
Subcritical Hopf bifurcations yield more complex post-bifurcation dynamics in the 
PDE wavetrains. In special cases where the subcritical bifurcation marks the end of 
the regime of stability, the post-bifurcation behavior in the spatial ODEs is chaotic, 
corresponding to wavetrains of the original PDEs which are spatially coherent, but have 
chaotic temporal dynamics. All the models are integrated numerically to investigate 
the post-bifurcation dynamics and chaotic regimes are characterized by computing 
power spectra, autocorrelation functions, and fractal dimensions. 



1 Introduction 

Morphogenesis or the occurrence of spatial form and pattern evolving from a spatially ho- 
mogenous state is a fundamental problem in developmental biology. A seminal contribution 
to this problem was made by Turing [1] who studied reaction-diffusion equations of the form 
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In [1] , the reaction functions (or kinematic terms) Ri and R2 were polynomials. However, the 
fundamental, and somewhat surprising, result that diffusion could destabilize an otherwise 
stable equilibrium leading to nonuniform spatial patterns (referred to as prepattern) is not 
dependent on particular forms of Ri and R2. 

The Turing instability in reaction-diffusion models thus provided a plausible and robust 
mechanism for the establishment of spatial prepattern, which could then generate biolog- 
ical patterns for gene activation. Numerous extensions and applications followed. These 
include early theoretical and analytical extensions [1-3]. In particular, Segel and Jackson 
[4] showed that spatial patterns may occur via Turing instability in macroscopic (extended 
Lotka-Volterra) models in population biology as well, particularly for species dispersing at 
different rates. They also provided a lucid physical explanation of how diffusion could in- 
deed generate instability, contrary to its usual interpretation as a smoothing mechanism. 
Applications in development biology were stimulated by the work of Meinhards and Gierer 
[5-7], primarily consisting of numerical simulations of reaction-diffusion systems in vari- 
ous geometries. Analytical work has confirmed and extended the results of [5~7], including 
bifurcation analysis and investigations of nonstationary (traveling-wave) patterns, spirals, 
solitary peaks, and fronts [8-12]. These are reviewed in [13]. Other work has focused on 
explaining the properties of spatial patterns [14-17] on the basis of chemical interactions 
and geometric considerations. Alternative explanations of pattern-formation, not based on 
reaction-diffusion equations and the Turing mechanism, have also been investigated [18]. 
Recent reviews of these and other related work on spatial pattern formation are given by 
Levin and Segel [19], Murray [20] and Edelstein-Keshet [21]. 

In order to incorporate various realistic physical effects which may cause at least one of 
the physical variables to depend on the past history of the system, it is often necessary to 
introduce time-delays into the governing equations. Factors that introduce time lags may 
include age structure of the population (influencing the birth and death rates), maturation 
periods (thresholds), feeding times and hunger coefficients in predator-prey interactions, 
reaction times, food storage times, and resource generation times. Models incorporating 
time delays in diverse spatially-homogenous biological systems are extensively reviewed by 
MacDonald [22], and in the context of predator-prey models, by Gushing [23]. These include 
continuous models such as the Kolmogorov, May, HoUing, Hsu, Leslie, and Caperon models, 
as well as discrete models. 

Consider (1.1) for the general two-species predator-prey model [24] with 

R^{N,P) = NF(N) -aNP — (1.2) 

k 

R2{N,P)^-PG{P)+^NP, 

where N{t) and P{t) are the prey and predator populations, respectively, e is the birth rate 
of the prey. A; > is the carrying capacity, a is the rate of predation per predator, and (3 is 
the rate of the prey's contribution to predator growth. 

In this paper, we initiate a fresh and detailed investigation of traveling spatial wave 
patterns of (1.1). In particular, we shall investigate in detail wavetrains with periodic and 
chaotic spatial variation. Toward this end, we consider traveling wave solutions of (1.1) in 
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the form of 



N{x,t) = N{0 (1.3) 
P{x,t) = P{C), 

where = x — vt is the travehng wave, or "spatial", variable, and v is the translation or 
wave speed, which will act as our bifurcation parameter. Substitution of Eqns. (1.2), (1.3) 
in (1.1) leads, after some simplification, to the four-mode dynamical system 



N^M (1.4) 

1 / eN'^\ 
M ^ —i^-vM - NF{N) + aNP + — j 

P = Q 



Q=^^[-vQ + PG{P)-/3NP)y 



where the overdot denotes ^. 

Here, however, we will follow [24,28] to consider the stability of the equilibria and the 
Hopf bifurcations of (1.4) for general functions F{N) and G{P). This is done in sections §2 
and §3. In section §4 we consider (1.4) for specific choices of F{N) and G{P) to determine 
the regions of phase-space where the system is volume contracting (dissipative), or volume 
expanding (dilatory). Also, note that the function F(N) incorporates the prey birth rate, 
and similarly, the function G{P) incorporates the predator death rate and is chosen so that 
this rate increases with predator density P. Section §5 considers the stability of physically 
relevant equilibria and Hopf bifurcation points for specific parameter values and choices 
of F{N), and G{P). Possible chaotic regimes are also dehneated there. The systems are 
numerically integrated and chaotic regimes are characterized by computing power spectra, 
correlation function and fractal dimensions [12]. Section §6 summarizes the results and 
presents the conclusions. 



2 Linear stability analysis 

The equilibrium, critical or fixed points of the system (1.4) (only nontrivial points are relevant 
since both the predator or prey population can not be zero) are 

(iV„,M„,P„,Q„)=(^,0,^:<^^2^,0) (2.1) 

In this section we will consider the Turing bifurcations in general predator-prey systems by 
considering the system (1.4) for general F{N) and G{P), which incorporate the prey birth 
rate and predator death rate. For numerical purposes in, the functions F{N) and G{P) are 
subsequently chosen to be 

A. F{N)^e, G(P)=7, e = 0, 

B. F{N)^e, G(P)=7, e = e. 
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C. F{N) = ko, G{P) ^d + cP,e^e, 



D. F{N) = ko{l + f ), G{P) = d + cP,e = e, 

E. F{N) = i±f , G{P) = 7(1 + kP'), e = 0. 

For the remainder of this paper we shall refer to these cases as System A-E. System B is a 
modified Lotka-Volterra two species model with diffusion, and 7 being the death rate of the 
predator. Notice that qualitative features of such models have been considered earlier, for 
instance for the Kolmogorov model without delay [2] and the May model with delay [1]. 

Following standard methods of phase-plane analysis the Jacobian matrix of (1.4) evalu- 
ated at the fixed point (A^o, -^0, -fo, Qo) is 



V 





aPo+^-FiNo)-NoF'{No) 



D2 



1 











aNp 


-PNo+G(Po)+PoG'{Po) 
D2 



\ 


1 



V 

'D2 



(2.2) 
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The eigenvalues of this matrix satisfy the characteristic equation 

g{X) = + biX^ + hX^ + 63A + 64 = 
where hi with i = 1, 4 are given by 



hi 
h2 
hs 

&4 



D1D2 

a{pkv'^ - eD2Go + kD2F^Go) + £)i(eG'o - pkF)G'Q 

a(3kDiD2 

v{ - eaGo + akF^^Go + {eGo - /3kF)G'o) 

al3kDiD2 
{l3kF^ - eGo) {af3k + (e - kFo)G'o)Go 
al3^kWiD2 



(2.3) 



(2.4) 



where Fq = F{No), Gq = G{Po), = ^^\No, and G'^ = ^^\Po- The Routh-Hurwitz 
criteria [20], giving the necessary and sufficient conditions Re{Xj) < 0, i = 1, 4 for stability 
of the steady state {Nq, Mq, Pq, Qq) are 



&i >0, 
64 > 0, 
- &3 > 0, 

6l(&2&3 - &I&4) - 63 > 0. 



(2.5a) 
(2.5b) 
(2.5c) 
(2.5d) 



Hence, instability of the steady state may arise for some traveling wave speed v if any of the 
above conditions are violated, i.e. 



hi<0, 64 <0, 6162 -&3<0, hi{h2h3 - hihi) - hi < 0. 



(2.6) 
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In other words, a regime of stability/instability could be created by varying the bifurcation 

parameter around the critical values v^. Note that the last condition (2.5d) corresponds, 
at equality, to a Hopf bifurcation with two roots of (2.3) having purely imaginary complex 
conjugate values. This condition is quartic in v and has the form 

f{v)^v'^{Av'^ + B)^Q (2.7) 

where, 

, _ (£>i + D2){- eaGo + akFl^Gp + (eCp - PkF)G',) 

B = ^,^,^^^3^3 {a'Go{e'DlGo + pk{D, + D,f{eG, - /3kFo)) 

+ {akD^GoF'^ + Di{l3kFo - eG)G'Q) ( - 2eaD2GoakD2G^F^ + L'i(/3A;Fo - eG)G'^) 



(2.9) 

The existence of real nonzero roots of (2.7), requires the necessary condition AB < 0, since 
on the Hopf curve 

% = Ty^. (2.10) 

The velocity of the traveling wave will give the change of stability of the steady state 
(A^O) -^0) Pq-i Qo) in the following manner: 

(a) If A > 0, and B < then the fixed point is stable in the region v e (— oo, V-) U {v+, oo) 
and unstable for v e {v-,v+), 

(b) if ^4 < 0, and -B > then the fixed point is unstable in the region v e (— oo, v_)U(v+, oo) 
and stable for v e (v-, v+), 

(c) if ^4 > 0, and B > then f{v) > 0, the steady state is stable Vv, 

(d) if A < 0, and B < then f{v) < 0, the steady state will be unstable \/v. 



3 Hopf bifurcation analysis 

We will perform a Hopf bifurcation analysis [4,6,9, 10] to show that as the value of v passes 
through the critical values v^, periodic solutions occur. To determine the behavior of the 
eigenvalues A as v varies , we use the following lemma. 

Lemma 1. The characteristic equation (2.3) with hi e 3?, n roots for i — 1,...,4 and 
discriminant given by 

A = n (r. - r,)^ (3.1) 
has a pair of purely imaginary roots and two real roots only when A < 0. 
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Proof. At V — v^:, b4 — ^ — ^, then the discriminant becomes 

Mbl - + 463)(&^&i + bib, - 4616263 + 46i)^ 
A- . (d.^j 

Since 63 > then A < when 6^ — 46162 + 463 > 0. We wiU see next why this last condition 
is satisfied. 

Rewriting (2.3) as 

(3.3) 

v( - eaGo + afcF^Go + (e^o - /3kF)G'o) {/3kFo - eGo) (a/3 A; + (e - A;Fo)G'(j)G'o 
^ aPkDiD^ ^ ^ a/3^kWiD2 ' 

and evaluating this on the Hopf curve, i.e., at v — yields to 

^( A, v^) = ( + u') ( A^ + sX + p), (3.4) 

where, 

ri,2 = l^ico = =Fm/ (3.5) 



6^-46162 + 463- 



and bi = bi{vzf). Since wc require that 61 > by (2.5a), and r,^4 G $R, then 6^—46162+463 > 0. 
This leads to a Hopf bifurcation setting as evidenced by the pair of imaginary eigenvalues 
ri^2 that oscillate with frequency 



G'„{eGo-PkF„) + aGo(kFi-e) 
" = V amO^ + D,) ■ 

□ 

In order to introduce the relevant notation, we state the Hopf bifurcation theorem. 
Theorem 1. Let 

^ = F(f,/i) (3.7) 

he an autonomous system of differential equations for each value of the parameter /i G 
(— /xcA^o)? where /iq is a positive number and the vector function F e C^{D x (— /xo,Aio)); 
where D is a domain in 3?". Suppose that the system (3.7) has a critical point Xo{ii), that 
is, 

F(fo(/i),/i) =0. (3.8) 

Let J{ijl) he the Jacohian matrix of system (3.7) at xq{^). Suppose that det{J{n) — XI) = 
has a complex conjugate pair of solutions A(/i),A*(/i) such that for fi > 0, Re\{fi) > 0; 
/X = 0, Re\{n) = 0; n < 0, Re\{n) < 0; and ^^^^^^|^=o > 0. Assuming that all other 
X's are distinct, (3.7) has a periodic solution in some neighborhood of /j, = and x in some 
neighborhood of xq{ij). 
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Proof. In order to apply this theorem to (3.3), we define the bifurcation parameter 

, 3.9 

V Vq 

then 

v{pi) = (3.10) 

with fi = ai V = Vq. For > 0, < fo and Re\{fi) > 0. For = 0, = fo and Re\{fi) = 0, 
and for fi < 0, v > vq and ReX{fi) < 0. Thus, the first set of conditions in the theorem are 
vahd, and it remains only to show that 

For Vq = Vzf, then 5'(A(i'o), I'd) = 9{l^i^^,vo) = by (3.4). Hence, implicitly differentiating 
g{X{v),v)^0, (3.3), yields: 



dX _ g _ uj{b'^ - u% + iuh'^) 



dv g 2a;(2a;2 - 62) 3a;26i) ' 



(3.12) 



where h\ — ^|^=^g, therefore 



dReX{vo) o-^^f 



4a;2(2a;2 - 62)^ + (63 - Su^iY ' 

where 



(3.13) 



^ = 2(2^2 _ _ ^ ^/^^^^ _ (314) 

Evaluating the required derivatives of 6i's at and using (3.6) yields 



4^;g G'o(-eGo + ^feFp) + aGoj-kF;, + e) 



Thus, assuming that Di + D2 > 0, and using (3.6), then ^' < 0, and hence < Q. 

All the conditions of the Hopf bifurcation theorem are satisfied. Thus, Hopf bifurcations 
occur and periodic solutions will exist in the neighborhood of Vq. □ 



4 Contracting/Dilatory behavior 

The stability of the bifurcating closed orbits may be investigated for each the specific choice 
of nonlinearity F{N) and G{P), by reducing the system to the center manifold (since one 
has two purely imaginary eigenvalues) as done in [28]. This will not be considered in here. 
Instead, we shall consider numerical solutions of (1.4) in the following section, which will 
allow both the verification of the preceding analysis and also yield more quantitative results. 
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We will concentrate on the five specific choices of F{N) and G{P) referred to us in this 

paper as systems A-E. For all models, the local rate of change of volume of the (A^, M, P, Q) 
phase-space in the vicinity of the fixed points (Nq, Mq, Po,Qo), which gives the local loga- 
rithmic rate of change of {N, M, P, Q) phase-space volume V is given by the trace of the 
Jacobian matrix of (2.2) at the fixed points, where Tr{J) = = — = — A 

necessary condition for the stabihty of the steady state is that 61 > by (2.5a), therefore 
models that start from stable/unstable fixed points (depending upon one or more of (2.5b)- 
(2.5d) is violated) will be locally dissipative, i.e., (phase-space volumes contract), so we may 
anticipate that the orbits may go to an attractor at infinity if the dissipation is weak, or 
dilatory (volumes expand) if (2.5a) is violated. If the fixed point is stable, the predator 
population is ultimately decimated, i.e., k and the rate of conversion (3 of prey into predator 
are not large enough to sustain the predator population. If the fixed point is unstable, for 
a parameter regime where the system is strongly dissipative, one might anticipate possible 
bounded chaotic dynamics evolving on a strange attractor. This will be tested numerically 
in the next section. 



4.1 System A 

Using F{N) = e, G{P) = 7, and e = 0, (1.4) becomes 

N = M (4.1) 

M ^ ^{-vM -eN + aNP) 
P = Q 

Q = -^{-vQ + ^P-l3NP), 

with equilibrium points (A^o, M), -fo, Qo) — 0> '^^^ characteristic equation (2.3) has 

coefficients 



bo — 

67 



b. 



'4 



D1D2 

Therefore, the Hopf curve (2.7) is 

(Di + D2)'e7 2 _ 
fiv) = ^3^- = 0' (4-3) 

hence, the bifurcation parameter is 

= 0. (4.4) 

The characteristic polynomial (2.3) evaluated at the fixed point and on the Hopf curve (4.4) 
has the form 

g{X,v^)^X^ + h, (4.5) 

and 64 = 64 (fq:). 
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4.2 System B 

Using F{N) = e, G{P) = 7, and e = e, (1.4) becomes 

N = M (4.6) 

M= -^(-vM-€N + aNP+^) 
Di^ k ' 

P = Q 

Q^^{-vQ + ^P-PNP), 

with equilibrium points {Nq, Mq, Pq,Qq) = (^,0, ^(1 — ^),0). The characteristic equation 
(2.3) has coefficients 



b2 



pkDiD2 

e'yilSk - 7) 



f3kDiD2 

Therefore, on the Hopf curve, the bifurcation parameter is 



e^Dl + pkjDr + D^n^-Pk) 

pm+D2) ' ^'-'^ 



where 

,2 



^^feS^i^i ( - + ^2)^' + 67i^2' + imDi + D,f - p'k\D^ + D,f). (4.9) 



The characteristic polynomial (2.3) evaluated at the fixed point and on the Hopf curve (4.8) 
has the form 

^(A, v^) = + h^X" + + 63A + 64, (4.10) 

where hi — hi{v^). 
4.3 System C 

Using F{N) = A;o, G{P) ^d + cP, and e = e, (1.4) becomes 



TV = 


M 




1 


M = 








P = 


Q 




1 


Q = 





(4.11) 
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with equihbrium points (No, Mq, Pq, Qo) = (^SS^' 0' ' 0)- The characteristic equa- 

tion (2.3) has coefficients 



&2 
&4 



q;(^A;^;2 - edL>2) + c{edDi - {eD2 + pkDi)ko + ev^) 

{ec + apk)DiD2 
v{ed{c — a) — c(e + ^k)) 

{ec + a(5k)DiD2 
{ad + cko){f3kko — ed) 



{ec + aPk)DiD2 
Therefore, on the Hopf curve, the bifurcation parameter is 



Ic'^iedDi + ekoD2 - ^kkoDif + a^d{deWl + /3k{ed - /3kko){Di + Da)^) + 9 



(ec + aPk) {Dl + D2) {ed{a - c) + c(e + /3k)ko) 

(4.13) 

where 



f{v) = j^^^—^^^^^[{ec + a/3k){Di + D2){{aD + cko){ed-/3kko){Di + D2) (4.14) 

{ed{a — c) + c(e + /3k)ko)) {c{/3kko — ed)Di + e{ad + cko)D2 — (ec + a/3k)v'^) . 

ec + a/? A; 

- DiD2(ec?(a - c) + cA;o(e + /3A;))') , 

and e = ac(2e^d'^DiD2 + edko{2eDj + /^^(D^ ^ ^2^^ _ /52p^2p^ ^ ^^^2^ _ 

In this case, the characteristic polynomial (2.3) evaluated at the fixed point and on the 
Hopf curve (4.13) has the same form as (4.10). 

4.4 System D 

Using F{N) = ko{l + f ), G{P) = + cP, and e = e ,(1.4) becomes 

TV = M (4.15) 

1 f .... . .r , „,.rn , {^-ko)N'\ 



M = —( -vM - koN + aNP + 

Di\ k 

^(^-vQ + dP- pNP + cP'), 



with equilibrium points {No, Mo, Po, Qo) = ( ("igt+^fc , 0, f^|^g^, 0). This system is not 
quantitatively different form System C, therefore all the equations (4. 12), (4. 13), and (4.14) 
will stay the same as long as we replace e ^ e — ko- 
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4.5 System E 

Using, F{N) = G{P) = 7(1 + kP^), and e = 0, (1.4) becomes 



(4.16) 



M 



D 
P = Q 



Q = ^(^-vQ + -fP{l + kP^)-pNPy 



The equilibrium points are found numerically by solving the following system which involves 
a quintic algebraic equation in Nq. 



No^jil + kP^) 



7 

Mo = 

Qo = o. 



(4.17) 
(4.18) 



1 + SNo 



The characteristic equation (2.3) has coefficients 
(L>i + D2)v 



hi 
b2 



D1D2 



(4.19) 



D1D2 



+ 



D2 



(1 + N^Y 

V 



+ pNoDi - 7^1(1 + 3kP^) 
(1 - aPo - No{No + aNoPo{2 + N^) - 2S)) 

PN^ + (7 + Po{a + 3kjPo))N^ + + 1 - Po(« + 3A;7Po) 



' D^D2{l + NlY 
- 7 + No{P + 25) - Nl (1 + 27 + 2Po(« + 3A;7Po)) 



64 = 



1 



(/3iVo - 7(1 + SfcPo^)) (gPo - 1 + iVo(iVo -25 + aiVoPo(2 + iVp^))) 

(1 + N^,y 



Therefore, on the Hopf curve, the bifurcation parameter Vzs^ can only be found numerically 
by (2.10), where 



D1 + D2 



PN^ + (7 + Po(a + 3fc7Po)X + 2/3iV3 + 1 - Po(q; + 3A;7Po) (4.20) 



D\Dl{l + NiY 
7 + 7Vo(;5 + 25) - Nl{l + 27 + 2Po(q; + 3A;7Po)) 

i^3^i(l + iV2)^ 



Di{l + iVo) T^i + 2^1^2(1 + iVo)^2 + i^2^; 
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and 



Ti = P^Nl + (7 + ^ikP^f - /3iVo(27 + Po(a + 67A;Po)) (4.21) 
T2 = a7Po(l + 3A;Po')A^o - /^^o + (2/35 + (1 + 3A;p2)(l + 2aPo))iVo' 

+ {(3- 275(1 + 3A;p2))iVo - 7(1 + 3A;P2)(1 - aPo) 
T3 = -a(3PoN^ + a^P^K - ^oi^PoNl + 2«Po(l + 2aPo)N^ - 2«(3/3 + 26)PoN^ 

+ (1 + 2aPo(l + 3aPo))Aro4 - 4(5 + ^(^g + 26)Po)N^ + 2(25^ _ 1 _ «(! _ 2aPo)Po)Aro' 

+ (4(5 - aPo(/3 + 4(5)) A^o + (1 - aPoY- 

5 Numerical results 

For the numerical results we will concentrate on our five systems, choosing for each system 
specific parameters that will show the dissipative or dilatory behavior. 

5.1 System A 

We choose parameters such that 64 = < 0. Since (2.5b) is violated it means that we 
start from an unstable fixed point in a constant volume space. In this case, f{v) > 0,Vf, 
therefore the steady state remains stable since both populations will annihilate. For the 
system parameters of = 1.5, e = 1, (3 = —2.75, 7 = —1.5, Di = 1.25, D2 = 2.1, and the 
bifurcation parameter f = 0.1, then 64 = = —0.571429 < 0, and the populations start 
to oscillate from the stationary point (A^'o, Mq, Pq, Qq) = (0.55, 0, 0.67, 0) with frequency 
CO = a/— 64 = 0.755. Because the space is contracting, since 63 = 0.128, eventually both 
predators/prey populations will assimilate each other and reach the stable equilibrium null 
populations. This attenuating behavior is presented in Fig.l. Note the stable periodic 
oscillations on the stable limit cycle created by a supercritical Hopf bifurcation at v = v^p = 0. 
If we were to increase v then the population would have terminated much faster. 




20 40 60 80 100 



Figure 1: Periodic time series for populations of System A, v = 0.1 
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5.2 System B 

For the system parameters given by the set a = —1.2, e = —3, /3 = —2, 7 = —2, Di = 1, 
D2 = 2, k = 2, the Hopf velocity is f=p = =f2. Therefore, the populations start to oscillate 

from any stationary point with frequency u = = To find the regimes when the fixed 

point changes stability, we find the coefficients of the Hopf condition (2.7), which are A = j^, 
and B = — |, and we analyze f{v). Hence, the fixed point {Nq, Mq, Po,Qo) = (1,0,1.25,0) 
is stable in the region v G {—00, —2) U (2, 00) and unstable for v G (—2, 2). 

Since the volume of the system is expansive on = —2, and contractive on f+ = 2, 
as we vary v around f=p we will expect different behavior on both sides of the bifurcation 
parameter. In a contracting space, v+ = 2, then bi = 3, 62 = 3.5, 63 = 1.5, 64 = 1.5, and hence 
the population will oscillate from any fixed point toward the equilibrium {Nq, Mq, Pq, Qo) = 
(1,0,1.25,0). 

When V = 2.2, the fixed point remains stable, hence the populations dissipate as in 
case A, but instead of reaching the null populations they will converge towards nonzero 
equilibrium values. This behavior is shown in Fig. 2. 

If f = 1.9, the fixed point becomes unstable, and, after an initial transient, both popu- 
lations settle onto the stable limit cycle created by the supercritical Hopf bifurcation. The 
corresponding spatially periodic wavetrain in spatial variable ( is shown in Fig. 3. 



N[e],p[S] 




50 100 150 200 



Figure 2: Periodic evolution for populations of System B, v = 2.2 

By contrast, on the left side of the bifurcation parameter, the system is expansive or 
dilatory at f_ = —2 and undergoes a subcritical Hopf bifurcation which occurs at f = 
f_. This corresponds to an unstable periodic orbit coexisting with an unstable fixed point 
(A^o, ^0, -Pq, Qo) = (1; 0, 1.25, 0), since (2.5c) is violated. For this case, bi = —3, 62 = 3.5, 
63 = —1.5, 64 = 1.5. Because the system is expanding then the only possibility is to have an 
attractor at infinity. Hence, the populations blow at a finite value of (. 



13 



N[e],P[e] 




Figure 3: Periodic evolution for populations of System B, v = 1.9 

5.3 System C/System D 

Since these two systems are similar as explained in previous section, for numerical simulations 
we will describe the behavior of only System D. Choosing the system parameters given by 
the set a = 1.25, e = 1, /3 = 2, c = 0.5, d = 2, k = 2, ko = 2, Di = 1 and D2 = -2, the 
equilibrium point is {Nq, Mq, Pq,Qq) = (1.55,0,2.22,0), while the bifurcation parameter on 
the Hopf curve is = ^5.03. Here, h = -2.51, 62 = -11.22, 63 = -0.83, 64 = -3.88. 
For these values, (2.5c) is not violated but (2.5a) and (2.5b) are, hence the fixed point is 
unstable. To find the regimes when the fixed point changes stability we find the coefficients 
of the Hopf condition (2.7), A = —0.04, and B = 1.05, and we analyze f{v). Hence, the 
fixed point will remain unstable in the region v G {—00, —5.03) U (5.03, 00) and stable for 
V G (—5.03,5.03). Within the stable region if f = —0.2, the volume is weakly expanding 
hence we anticipate that the obits may go to an attractor at infinity. From Fig. 5 we can 
see the aperiodic behavior of the populations. The orbits fly off to an attractor at infinity 
as shown in Fig. 4 by both N{Q, and -P(C) blowing up around ( = 96. 

If t> = 1.2 the volume is dissipative but the fixed point is unstable, hence the populations 
experience qusiperiodic behavior or bounded chaotic behavior. We will present this case 
next. 

5.4 System E 

Since in this case the fixed points can not be found analytically, due to a quintic algebraic 
equation, we will solve this case completely numerically. For the parameters set a = 1.7, 
(3 = —2.1, 7 = —2, 6 = 0.6, k = —2, Di = —1 and D2 = 2, the equilibrium point is 
{Nq, Mq, Pq,Qq) = (0.19,0,0.63,0), while the bifurcation parameter on the Hopf curve is 
= ^1.8. Here, 61 = -0.89, 62 = -3.252, 63 = 2.84, 64 = 0.27. For these values, (2.5c) is 
not violated but (2.5a) is, therefore the fixed point is unstable. To find the regimes when the 
fixed point changes stability we find the coefficients of the Hopf condition (2.7), A = 0.39, 
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Figure 4: Aperiodic evolution for populations of System D, f = —0.2 




Figure 5: Attractor at infinity for populations of System T), v = —0.2 

and B = —1.027, and we again analyze f{v). Hence, the fixed point will become stable 
in the region v G (— oo,— 1.8) U (1.8, oo) and stable for v G (—1.8,1.8). For v = —1.1 the 
volume is dissipative and as explained above the populations will behave chaotically. Fig 
6 shows the numerical solutions for A^(C) and -P(C) vs. the spatial variable (. Notice the 
strange aperiodic dynamics. Note that unlike Fig. 4 the solution remains unbounded for 
all time. The 3D phase space plot in the (A^, M, P) space is shown in Fig. 7. Notice that 
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the solutions retrace the same region of phase space repeatedly, suggesting bounded chaotic 
dynamics on an attr actor. 




Figure 7: Attractor for populations of System E, v = —1.1 

In order to confirm this and further characterize the suspected chaotic solutions, we 
employ the standard numerical diagnostics [25,26] i.e., the power spectral density, the au- 
tocorrelation function, and the fractal dimensions. The power spectral density and the 
autocorrelation function of A^(C) are computed using codes from "Numerical Recipes in C" 
[27], and the former is shown in Figs. 8, 9. The "broad" peaks in the power spectral density 
plot are indicative of chaos and randomness. 
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Figure 8: Power spectral density of System E, v = —l.l 
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Figure 9: Log PSD vs. frequency of System E, v = —1.1 

However, we move on to a more quantitative and definitive numerical diagnostic, i.e., 
the fractal dimension [28]. In order to distinguish low-dimensional (deterministic) chaos 
from strong randomness, one computes the dimensions as discussed below. Of several pos- 
sible alternative definitions [25, 26] for the fractal dimensions, we employ the cluster fractal 
dimension D of Termonia and Alexandrowicz which is defined by 

n = fc[i?(n)]^,n ^ oo, (5.1) 

where R(n) is the average radius of an E-dimensional ball containing n points. Thus, D is the 
slope of the plot of logn vs. logR(n). More usefully, if a scaling law (5.1) exists it would show 
up as a horizontal line on a plot of dlogn/ dlogR{n) vs. logn with the height of the line being 
a measure of D. Fig 10 shows D which is the height of the approximate horizontal straight 
line, and we may estimate the converged cluster fractal dimension to be approximately 1.6. 
This confirms that the System E indeed possesses bounded low dimensional (deterministic) 
chaotic solutions evolving on a strange attr actor with dimension D. 
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Figure 10: Cluster dimensions calculation of System E, v = —1.1 

6 Summary and conclusions 

In this paper traveling wave pattern formation in general reaction-diffusion/predator-prey 
models including diffusion in the interspecies interaction terms has been considered. For 
our first two specific choices of nonlinear terms, the numerical and mathematical results 
presented here show either stable equilibrium behavior as in System A, or stable periodic 
spatial patterns as in System B. Systems C/D exhibit aperiodic spatial behavior (including 
a finite-time singularity using ODE terminology, or an attractor at infinite in dynamical 
systems parlance). For System E we also have aperiodic behavior within a diffusive volume, 
hence the patterns evolve chaotically on a strange attractor. 

Various immediate applications of these results suggest themselves. In particular, fu- 
ture work will address specific reaction-diffusion systems such as the Belousov-Zhabotinsky 
system. Other work in progress includes pulse-train dynamics, as well as the possibility of 
unsteady pulse solutions in such systems, similar to those recently observed and analyzed in 
the famous cubic-quintic Ginzburg-Landau equation. 
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